from vectors import *
from visuals import *
from PointEntity import *
from VoxelMap import *

from Tkinter import *
import numpy
from numpy.random import poisson
from random import randrange
import simplejson as sj
import sys

class Simulation:
	"""The class which binds the heatmap, objects, and entities together in
	a single system.
	"""
	next_floor_id = 0
	next_tile_id = 0
	next_detector_id = 0
	next_emitter_id = 0
	next_obstacle_id = 0
	next_node_id = 0
	
	time = 0.0
	window_geometry = "400x400"
	canvas_dim = vec2(400, 400).int()
	timestep = .1
	visual_timestep = 1.0
	time_until_visuals_update = 0.0
	duration = float("inf")
	source_detected = False
	simulation_running = True
	assumed_background = 8
	
	visuals = None
	heatmap = None

	# Obstacle Management
	obstacles = {}
	obstaclemap = None
	obstacleLocators = {} # Keys are by Detector ID, Assuming single Source
	
	floors = {}
	detectors = {}
	detectors_by_hits = []
	max_group_size = 4
	max_group_distances = {}
	# key:value :: group_size:[ID_tuple1, ID_tuple2...]
	groups = {}
	groups_hits = {}
	hottest_groups = {}
	hottest_groups_hits = {}
	
	emitters = {}
	nodes = {}
	node_graph = {}

	def __init__(self, \
	             map_file, \
	             debug_mode=True, \
	             heatmapOn=True, \
	             visualsOn=True, \
	             outputfile=None,
	             visualsfile="default"):
		""" """
		temp = str(self.__class__)
		self.class_string = temp[temp.find('.')+1:]
		MSG = self.msg
		self.debug_mode = debug_mode
		MSG("Debug mode is active!  Expect a wall of text." + "\n")#
		MSG("Conducting a %s" % self.class_string + "\n")#
		MSG("Setting visual mode to %s" % visualsOn + "\n")#
		self.map_file = map_file
		self.heatmapOn = heatmapOn
		# can't have visuals without a heatmap to drive it
		self.visualsOn = visualsOn and heatmapOn
		self.outputfile = outputfile
		self.outputfile_contents = ""
		self.visualsfile = visualsfile
		self.load(map_file)
		if self.visualsOn:
			self.visuals = Visuals(self)
		if self.outputfile:
			if self.class_string == "KSigmaSimulation":
				self.outputfile_line = "%f\t" * (self.max_group_size + 1) + "\n"
			elif self.class_string == "SweepSimulation":
				self.outputfile_line = "%f\t%f\t%f" + "\n"
			else:
				raise RuntimeError("something has gone badly")

	def msg(self, message):
		"""A message that will only be printed if debug_mode is set to True."""
		if self.debug_mode: sys.stderr.write(message)

	def load(self, map_file):
		"""Given a string representing a map file, attempts to load the state
		of the environment to this simulation.
		map_file is a text file specifying the environment in which the
		test will be conducted.  It has a specific format which makes use of
		the JSON method of encoding, and should be evident from the example
		map files."""
		MSG = self.msg
		try:
			MSG("Attempting to load map file: " + map_file + "\n")#
			f = open(map_file, "r")
		except IOError:
			raise IOError("Couldn't access map file: " + map_file)
		contents_list = f.readlines()
		contents_string = ""
		for line in contents_list:
			if "#" in line:
				line = line[:line.find("#")] + "\n"
			contents_string += line
		MSG("Attempting to extract data from JSON map file" + "\n")#
		mapdata = sj.loads(contents_string)
		
		MSG("Attempting to load general information" + "\n")#
		self.window_geometry = "%sx%s" % tuple(mapdata["GENERAL"]["displaywindow"])
		self.canvas_dim      = vec2(mapdata["GENERAL"]["displaywindow"])
		self.timestep        = mapdata["GENERAL"]["timestep"]
		self.visual_timestep = mapdata["GENERAL"]["visual_timestep"]
		self.duration        = mapdata["GENERAL"]["duration"]
		self.max_group_size      = mapdata["GENERAL"]["max_group_size"]
		MSG("Processing maximum group distances")
		self.load_max_group_distances(mapdata)
		self.background = mapdata["GENERAL"]["background"]
		self.assumed_background = mapdata["GENERAL"]["background"]
		
		if self.heatmapOn:
			MSG("Attempting to create a HeatMap" + "\n")#
			if self.class_string == "KSigmaSimulation":
				self.heatmap = KSigmaHeatMap(mapdata["HEATMAP"]["corner1"],
											 mapdata["HEATMAP"]["corner2"],
											 mapdata["HEATMAP"]["voxelsize"])
			elif self.class_string == "SweepSimulation":
				self.heatmap = BumpSweepHeatMap(mapdata["HEATMAP"]["corner1"],
											mapdata["HEATMAP"]["corner2"],
											mapdata["HEATMAP"]["voxelsize"])
			else:
				raise RuntimeError("No heatmap type is associated with a " + \
				                   "simulation of type: " + class_string)
			self.heatmap.sim = self
			self.heatmap.initialize_likelihood_array()
			MSG("%dx%dx%d heatmap initialized." % \
				tuple(self.heatmap.resolution) + "\n")#
		else:
			MSG("HeatMap is turned off.")

		if self.visualsOn:
			MSG("Visuals are on: attempting to load data from FLOORS." + "\n")#
			for floor in mapdata["FLOORS"]:
				self.add_floor(floor)
			MSG("Defined %d floors in the environment." % len(self.floors) + "\n")#
		
		MSG("Attempting to load data from DETECTORS." + "\n")#
		for detector in mapdata["DETECTORS"]:
			self.add_detector(detector)
		MSG("Added %d detectors to the environment." % len(self.detectors) + "\n")#
		if self.heatmapOn:
			self.heatmap.detectors = self.detectors
			MSG("and linked them to the HeatMap" + "\n")#
		self.build_detector_groups(self.max_group_size)
		MSG("and built the list of detector groups." + "\n")#
		
		MSG("Attempting to load data from EMITTERS." + "\n")#
		for emitter in mapdata["EMITTERS"]:
			self.add_emitter(emitter)
		MSG("Added %d emitters to the environment." % len(self.emitters) + "\n")#
		
		MSG("Attempting to load data from OBSTACLES." + "\n")#
		for obstacle in mapdata["OBSTACLES"]:
			self.add_obstacle(obstacle)
		MSG("Added %d obstacles to the environment." % len(self.obstacles) + "\n")#
		
		MSG("Attempting to create a Obstacle Map" + "\n")#
		self.obstaclemap = ObstacleMap(mapdata["OBSTACLEMAP"]["corner1"],
		                       	mapdata["OBSTACLEMAP"]["corner2"],
		                       	mapdata["OBSTACLEMAP"]["voxelsize"],
								self.obstacles)
		MSG("%dx%dx%d obstacle map initialized." % \
				tuple(self.obstaclemap.resolution) + "\n")#
	
		self.build_obstacleLocator()
		MSG("Obstacles are located" + "\n")#
		
		MSG("Attempting to load data from NODES." + "\n")#
		for node in mapdata["NODES"]:
			self.add_node(node)
		MSG("Added %d nodes to the environment." % len(self.obstacles) + "\n")#

		MSG("Attempting to build nodegraph" + "\n")#
		self.build_nodegraph()
		
		MSG("All map data seems to be in order!" + "\n")#
		self.mapdata = mapdata

	def load_max_group_distances(self, mapdata):
		"""Loads the dictionary of group distances, making sure that its keys
		are integers and its values are floating point numbers."""
		try:
			temp = mapdata["GENERAL"]["max_group_distances"]
		except KeyError:
			MSG("Warning: No max group distances specified.")
		# cast as int and float
		for key, value in temp.items():
			try:
				self.max_group_distances[int(key)] = float(value)
			except:
				self.msg("Warning: bad pair encountered and ignored: " + \
				         "size %s, distance %s" % (key, value))
		# make sure there's no limit on unspecified group sizes
		for i in range(1, self.max_group_size + 1):
			if i not in self.max_group_distances:
				self.max_group_distances[i] = float('inf')

	def add_floor(self, floor):
		"""Adds a "floor" to the display.  A floor is a portion of the heatmap
		compressed to fit in two dimensions.  Its only real purpose is to provide
		the visualization of the probabilities in a volume.
		Floors are defined in terms of the slices of the heatmap they occupy, 
		because defining them in terms of volume is a pain in the ass.
		"""
		if not self.heatmap:
			raise RuntimeError("No heatmap specified yet, cannot define floors!")
		id = self.next_floor_id
		self.next_floor_id += 1

		try:
			self.floors[id] = Floor(id, \
			                        self.heatmap, \
			                        floor["nlo"][0], \
			                        floor["nhi"][0], \
			                        floor["nlo"][1], \
			                        floor["nhi"][1], \
			                        floor["nlo"][2], \
			                        floor["nhi"][2], \
			                        floor["scalemin"], \
			                        floor["scalemax"])
		except KeyError:
			print "Couldn't add floor %d, critical information was" + \
			      " missing from the map file." % id

	def add_detector(self, detector):
		"""Adds a detector to the set of detectors in this simulation.
		Note that this is NOT passed a detector, but a dictionary "shaped like" a
		detector.
		Primarily for use in the load method.  As such, it relies heavily on the
		existence of many fields in the detector's dictionary.
		Note that it relies on the next_detector_id data element to assign an id.
		"""
		id = self.next_detector_id
		self.next_detector_id += 1
		try:
			self.detectors[id] = Detector(id,
			                              detector["position"],
			                              detector["velocity"])
			d = self.detectors[id]
		except KeyError:
			print "Couldn't add detector %d, critical information was" + \
			      " missing from the map file." % id
			return
		if "name" in detector:
			d.name = detector["name"]
		if "size" in detector:
			d.size = detector["size"]


########## oh god fix this
		d.num_wedges = 360
		d.initializeWedgeDecays()
		d.setRotationalVelocity(10)
#		if ("wedges" in detector):
#			d.num_wedges = detector["wedges"]
#			d.initializeWedgeDecays()
#			if ("rotational_velocity" in detector):
#				if ("rotational_velocity_mode" in detector):
#					d.setRotationalVelocity(detector["rotational_velocity"], detector["rotational_velocity_mode"])
#				else:
#					d.setRotationalVelocity(detector["rotational_velocity"])

		self.detectors[id].sim = self
	
	def add_emitter(self, emitter):
		"""Adds a emitter to the set of emitters in this simulation.
		Note that this is NOT passed an emitter, but a dictionary "shaped like" i
		an emitter.
		Primarily for use in the load method.  As such, it relies heavily on the
		existence of many fields in the emitter's dictionary.
		Note that it relies on the next_emitter_id data element to assign an id.
		"""
		id = self.next_emitter_id
		self.next_emitter_id += 1
		try:
			position = emitter["position"]
			# Make sure that emitter is at least 1 meter away from any detector
			if len(position) == 2:
				c1 = vec3(position[0])
				c2 = vec3(position[1])
				position = c1.min(c2) + (abs(c1-c2)).random()
				while any(map( \
					lambda d: (self.detectors[d].p - position).mag() < 1,
					self.detectors)):
					position = c1.min(c2) + (abs(c1-c2)).random()
				
			self.emitters[id] = Emitter(id,position,emitter["velocity"])
		except KeyError:
			print "Couldn't add emitter %d, critical information was " + \
			      "missing from the map file." % id
			return

		
		if "name" in emitter:
			self.emitters[id].name = emitter["name"]
		# Strength in Millicuries. Defaults to 1
		if "intensity" in emitter:
			self.emitters[id].intensity = emitter["intensity"]
		self.emitters[id].sim = self
	
	def add_obstacle(self, obstacle):
		"""Adds the obstacles that are provided by the map file.
		Primarily for use in the load method."""
		id = self.next_obstacle_id
		self.next_obstacle_id += 1
		try:
			self.obstacles[id] = Obstacle(id, 
											obstacle["pt1"],
											obstacle["pt2"],
											obstacle["velocity"])
		except KeyError:
			print	"Could not add obstacle %d, critical information was" + \
					" missing from the map file." % id
			return
		if "name" in obstacle:
			self.obstacles[id].name = obstacle["name"]
		if "type" in obstacle:
			self.obstacles[id].type = obstacle["type"]
		else:
			print "Warning: obstacle %s has no shielding type, default to air." \
					% id

	
	def add_node(self, node):
		"""Adds a node to the set of nodes in the given simulation.
		Primarily for use in the load method."""
		id = self.next_node_id
		self.next_node_id += 1
		try:
			self.nodes[id] = GraphNode(id, 
										node["position"],
										node["velocity"])
		except KeyError:
			print 	"Couldn't add node %d, critical information was" + \
					" missing from the map file." % id
			return
		if "name" in node:
			self.nodes[id].name = node["name"]
	
	def build_obstacleLocator(self):
		"""Builds the set of all Obstacle Locators between every pair of
		emitter, detector pair. Must be called after the obstacle mapping has
		been completed. Locators are stored by Detector as keys"""
		for E in self.emitters:
			for D in self.detectors:
				self.obstacleLocators[self.detectors[D]] = \
						self.obstaclemap.getObstaclesBetPts(	\
											self.emitters[E].p, 
											self.detectors[D].p)
	
	def build_detector_groups(self, n=4):
		"""Builds the list of all single detectors, pairs, up to n-tuplets of
		detectors.
		Note: first clears the dictionary of groups.
		Note: also builds groups_hits as it goes."""
		if n > len(self.detectors):
			print ("Warning: maximum requested group size (%d) exceeds " % n) + \
			      ("of detectors (%d).  Lowering group size accordingly." % \
			      len(self.detectors))
			n = len(self.detectors)
		self.groups = {}
		for i in range(1, n+1):
			self.groups[i] = []
			self.groups_hits[i] = []
		for i in range(len(self.detectors)):
			self.build_detector_groups_helper([i], n)
		
	def build_detector_groups_helper(self, group, max_group_size):
		"""The recursive helper function that allows for construction of groups
		of arbitrary size."""
		detectors = self.detectors
		self.groups[len(group)].append(tuple(group))
		self.groups_hits[len(group)].append(0)
		if len(group) < max_group_size:
			max_dist = self.max_group_distances[len(group) + 1]
			# This means "for the rest of the detector IDs that haven't been
			# inserted into a tuple"
			for i in range(group[-1]+1, len(self.detectors)):			
				it_is_too_far_away = False
				if max_dist < float('inf'):
					# See if any of the detectors are too far from this new
					# detector.  Adding a detector too far away would create
					# an invalid group.
					for j in group:
						if detectors[i].distanceTo(detectors[j]) > max_dist:
							it_is_too_far_away = True
							break
				if it_is_too_far_away:
					continue
				else:
					self.build_detector_groups_helper(group + [i], \
				                                  max_group_size)

	def build_nodegraph(self):
		"""Builds the adjacency matrix after all of the nodes have been added"""
		
		# Initialized to be complete graph
		for node in self.nodes:
			temp = range(self.next_node_id)
			temp.remove(node)

		## TODO: Need a function to connect if no obstacles are present when
		## nodes are raytraced

	def update_detectors_by_hits(self):
		"""Makes the detectors_by_hits list of pairs up-to-date.  This is a
		list of (detectorID, hits) sorted in descending order."""
		self.detectors_by_hits = \
			[(self.detectors[id].hits, self.detectors[id]) \
			 for id in range(len(self.detectors))]
		self.detectors_by_hits.sort()
		self.detectors_by_hits.reverse()
		


class KSigmaSimulation(Simulation):
	
	k = {}
	kfusion = 0
	kthreshold = 0
	
	def update_groups_hits(self):
		"""Runs through each group of detectors and updates the groups_hits
		dictionary's lists with the total for each group.  Also updates
		the hottest_groups dictionary"""
		for n in range(1, self.max_group_size + 1):
			for i in range(len(self.groups[n])):
				# "the i'th group of size n"
				self.groups_hits[n][i] = 0
				for d in self.groups[n][i]:
					self.groups_hits[n][i] += self.detectors[d].hits
			max_hits = max(self.groups_hits[n])
			self.hottest_groups[n] = \
				self.groups[n][self.groups_hits[n].index(max_hits)]
			self.hottest_groups_hits[n] = max_hits

	def update_k_values(self):
		"""Updates k1, k2, k4, kfusion according to current group counts."""
		d_size = 6
		# New Poisson number might vary but for faster computation
		total_background = poisson(self.background * self.time * d_size)
		f1 = 1.01
		f2 = 1.1
		f4 = 1.5
		for i in range(1, self.max_group_size + 1):
			self.k[i] = float(self.hottest_groups_hits[i] - \
			                  i * total_background) \
		                     / (i * total_background) ** .5
			if self.k[i] < 0:
				self.k[i] = 0
			#print self.k[i]
			#print self.hottest_groups_hits[i]
			#print self.time
		# Change to OR of all Three 
		self.kfusion = self.k[1]/f1 + self.k[2]/f2 + self.k[4]/f4


	def is_there_a_source(self):
		"""If kfusion > kthreshold, we think there is a source.  Returns
		and sets source_detected to True or False accordingly."""
		self.source_detected = (self.kfusion > self.kthreshold)
		if self.source_detected:
			self.msg("Radioactive anomaly detected.")
		else:
			self.msg("No anomalies detected.")
		return self.source_detected

	def detect_emitters(self, timestep=1.0):
		"""Update each detector's hit count with the counts it should be
		getting from each emitter."""
		for detector_id in range(len(self.detectors)):
			d = self.detectors[detector_id]
			#	Debug
			#	print self.obstacleLocators[d]
			for emitter_id in range(len(self.emitters)):
				e = self.emitters[emitter_id]
				# d.hits += d.getRadCounts(timestep,e)
				# Obstacles in Place
				d.hits += d.getRadCounts(timestep, e)
				
				# print str(d) + ": " + str(d.hits)

	def detect_background(self, timestep=1.0):
		"""Given the mean value of background hits, updates the detectors with
		some number of hits centered around the mean"""
		for detector_id in range(len(self.detectors)):
			d = self.detectors[detector_id]
			hits = poisson(self.background * timestep * d.size)
			# Update the total hit values
			d.hits += hits

	def step(self, saveVisuals=False):
		"""Actually takes one step forward in time, updating all detectors'
		counts and seeing which group has the most hits, before executing the
		k-sigma algorithm.
		TODO: differentiate this from the timestep that would be relevant to
		motion."""
		self.time += self.timestep
		# first, get all the detectors up-to-date on radiation counts.
		self.detect_emitters(self.timestep)
		self.detect_background(self.timestep)
		self.update_detectors_by_hits()
		# also updates hottest_groups
		self.update_groups_hits()
		self.update_k_values()
		self.is_there_a_source()
		if self.heatmapOn:
			self.heatmap.update()
			est_pos = self.heatmap.estimatePosition()
			if self.outputfile:
				# Debugging
				# print est_pos, self.emitters[0].p
				# requires single source
				if self.emitters:
					dist_error = (est_pos - self.emitters[0].p).mag()
				else:
					dist_error = 0
				self.outputfile_contents += self.outputfile_line % \
			                            (self.k[1], self.k[2], self.k[3], 
										self.k[4],dist_error)
			if self.visualsOn:
				self.time_until_visuals_update -= self.timestep
				if self.time_until_visuals_update <= 0.0:
					self.time_until_visuals_update += self.visual_timestep
					self.visuals.update()
					if saveVisuals:
						filename = self.visualsfile + \
								   str(self.time).replace(".", "_") + \
								   ".ps"
						self.visuals.c.postscript(file=filename)
		if self.time >= self.duration:
			self.simulation_running = False

	def run(self, saveVisuals=False):
		"""Runs the simulation through to completion."""
		MSG = self.msg
		while self.simulation_running:
			self.step(saveVisuals)
			# Might want debug info here
		if self.visualsOn:
			self.visuals.tk.quit()
		if self.outputfile:
			f = open(self.outputfile, "w")
			f.write(self.outputfile_contents)
			f.close()



class SweepSimulation(Simulation):
	"""Makes use of directed Detectors.  Currently, does not have the detectors
	actually sweep out a region; instead, goes on the approximation of
	infinite sweep speed."""
	
	def detect_emitters(self, timestep=1.0):
		"""Updates the likelihoods of all positions in the heatmap."""
		hm = self.heatmap
		for d in range(len(self.detectors)):
			# we want the Detector, not its id.
			d = self.detectors[d]
			# note that vec3.angle normalizes this, so we don't
			# need to here.  to_p stands for "to point".
			for e in range(len(self.emitters)):
				# we want the Emitter, not its id.
				e = self.emitters[e]
				# to_e is "to emitter".
				to_e = e.p - d.p
				counts = d.getRadCounts(timestep * d.arc_fraction, e)
				d.hits += counts
				for i in range(hm.n.x):
					for j in range(hm.n.y):
						for k in range(hm.n.z):					
							to_p = hm.voxelPos(i,j,k) - d.p
							offset_angle = to_p.angle(to_e, d.mode)
							if offset_angle < (2 * d.theta):
								voxel_counts = counts \
									* (d.theta - offset_angle / 2) / d.theta
								hm.likelihood[i][j][k] += voxel_counts
	
#	def detect_background(self, timestep=1.0):
#		"""Given the mean value of background hits, updates the detectors with
#		some number of hits centered around the mean"""
#		wedges = self.heatmap.wedges
#		num_wedges = self.heatmap.num_wedges
#		for d_id in range(len(self.detectors)):
#			d = self.detectors[d_id]
#			# hope for degs!
#			wedges_in_arc = int(d.theta * 2 / (360 / num_wedges))
#			arc_mean = (self.background * timestep \
#			            * d.size * d.window_fraction) \
#			            / (num_wedges)
#			arr = numpy.empty(num_wedges + wedges_in_arc - 1)
#			for w in range(num_wedges):
#				arr[w] = poisson(arc_mean)
#			# create the illusion of wrapping around
#			arr[num_wedges:] = arr[:wedges_in_arc - 1]
#			total_background_this_step = 0#####
#			# for each wedge in this detector's list
#			for w in range(num_wedges):
#				# find the sum of poisson rolls
#				counts = sum(arr[w:w+wedges_in_arc])				
#				d.hits += counts
#				total_background_this_step += counts######
#				# for each voxel in the wedge
#				for (i,j,k) in wedges[d_id][w]:
#					# contribute those counts
#					self.heatmap.likelihood[i][j][k] += counts
#			print "%s saw %s background counts this round" \
#				% (d, total_background_this_step)
	
	def photonsHit(self, d, counts=1):
		"""
		Used to update the heatmap and the detector's hits count when n photons
		strike detector d in its current position.
		"""
		d.hits += counts
		hm = self.heatmap
		w_front = d.wedge
		w_back  = d.wedge - d.num_wedges / 2
		# first, add the full count to the wedge dead ahead.  (Yes, this
		# means the first element of wedge_decays, which should be 1, is not
		# useful.)
		# I'm concatenating the wedge ahead and the wedge behind.  In a twelve-
		# wedge system, if one were 12 o'clock, the other would be 6.
		for (i,j,k) in \
		  hm.wedges[d.id][w_front]:## + hm.wedges[d.id][w_back]:
			hm.likelihood[i][j][k] += counts
#		print "Wedges %3s, %3s: counts = %g" % (w_front, w_back % 360, counts)
		# in the aforementioned twelve-wedge system, the next four wedges to
		# concatenate would be 1, 11, 5, and 7 o'clock, and so on.
# #		for w in range(1, len(d.wedge_decays)):
# #			decayed_counts = counts * d.wedge_decays[w]
#			print "wedges %3s, %3s, %3s, %3s: decayed_counts = %g" % \
#				((w_front + w) % d.num_wedges, (w_front - w) % d.num_wedges, \
#				 (w_back + w) % d.num_wedges,  (w_back - w) % d.num_wedges, \
#				 decayed_counts)
# #			for (i,j,k) in \
# #			  hm.wedges[d.id][(w_front + w) % d.num_wedges] + \
# #			  hm.wedges[d.id][(w_front - w) % d.num_wedges]:# + \
##			  hm.wedges[d.id][(w_back  + w) % d.num_wedges] + \
##			  hm.wedges[d.id][(w_back  - w) % d.num_wedges]:
# #				hm.likelihood[i][j][k] += decayed_counts

	def detect_background(self, d, timestep=1.0):
		"""
		Strikes Detector d with background radiation probabilistically
		appropriate to its size, assuming current orientation.  Does not
		rotate the detector.  (use stepWedge)
		TODO: Verify.
		"""
		background_counts = poisson(self.background * d.size * timestep)
#		print "Background for wedge %3s: %i" % (d.wedge, background_counts) 
		self.photonsHit(d, background_counts)

	def detect_emitters(self, d, timestep=1.0):
		"""
		Strikes Detector d with radiation from each emitter on the field.
		Realize that regardless of the position of the emitter, the detector
		is set to believe that the most likely place for it to come from was
		normal to the surface.
		"""
		for e_id in range(len(self.emitters)):
			e = self.emitters[e_id]
			####### OBSTACLES DISABLED
			counts = d.getRadCounts(timestep, e, use_obs=False, use_cos=True)
			if d.id == 0:
				print "%s: " % d.wedge,
				print counts
			self.photonsHit(d, counts)
			

	def step(self, saveVisuals=False):
		"""Actually takes one step forward in time, updating all detectors'
		counts and seeing which group has the most hits, before executing the
		k-sigma algorithm.
		TODO: differentiate this from the timestep that would be relevant to
		motion."""
		self.time += self.timestep
		# SweepSimulation is super dependent on a heatmap.  Like, this won't
		# work without it.
		for d_id in range(len(self.detectors)):
			d = self.detectors[d_id]
			for substep in range(d.wedges_per_timestep):
#				print
#				print
#				print "BACKGROUND:"
	#			self.detect_background(d, self.timestep / d.wedges_per_timestep)
#				print "EMITTERS:"
				self.detect_emitters(d, self.timestep / d.wedges_per_timestep)
				d.stepWedge()

		# SweepHeatMap updates take care of the hits updating.
		est_pos = self.heatmap.estimatePosition()
		est_intensity = self.heatmap.estimateIntensity()
		if self.outputfile:
			# requires single source
			if self.emitters:
				dist_error = (est_pos - self.emitters[0].p).mag()
				intensity_error = est_intensity / self.emitters[0].intensity
			else:
				dist_error = 0
				intensity_error = est_intensity
			self.outputfile_contents += self.outputfile_line \
				% (dist_error, est_intensity, intensity_error)
		if self.visualsOn:
			self.time_until_visuals_update -= self.timestep
			if self.time_until_visuals_update <= 0.0:
				self.time_until_visuals_update += self.visual_timestep
				#subtractMean
				self.visuals.update(subtractMean=False)
				if saveVisuals:
					filename = self.visualsfile + \
					           str(self.time).replace(".", "_") + \
					           ".ps"
					self.visuals.c.postscript(file=filename)
		if self.time >= self.duration:
			self.simulation_running = False

	def run(self, saveVisuals=False):
		"""Runs the simulation through to completion."""
		MSG = self.msg
		self.heatmap.makeWedges(360)
		while self.simulation_running:
			self.step(saveVisuals)
			MSG("TIME: %g\n" % self.time)
			MSG("I think it's at %s\n" % self.heatmap.estimated_position)
			MSG("But it's at %s\n" % self.emitters[0].p)
			MSG("\n")
		if self.visualsOn:
			self.visuals.tk.quit()
		if self.outputfile:
			f = open(self.outputfile, "w")
			f.write(self.outputfile_contents)
			f.close()
